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ABSTRACT 

We present the simage software suite for the simulation of artificial extragalactic images, based empir- 
ically around real observations of the Hubble Ultra Deep Field (UDF). The simulations reproduce galaxies 
with realistic and complex morphologies via the modeling of UDF galaxies as shapelets. Images can be 
created in the B, V, i and z bands for both space- and ground-based telescopes and instruments. The sim- 
ulated images can be produced for any required field size, exposure time, Point Spread Function (PSF), 
telescope mirror size, pixel resolution, field star density, and a variety of detector noise sources. It has 
the capability to create images with both a pre-determined number of galaxies or one calibrated to the 
number counts of pre-existing data sets such as the HST COSMOS survey. In addition, simple options 
are included to add a known weak gravitational lensing (both shear and flexion) to the simulated images. 
The software is available in Interactive Data Language (IDL) and can be freely downloaded for scientific, 
developmental and teaching purposes. 

Subject headings: Simulations - Cosmology: weak lensing - Galaxies: Surveys 



1. Introduction 

In the next decade, the quantity of data available to 
cosmology will rapidly increase. New telescopes, both 
on the ground and in space, promise to image many 
thousands of square degrees. The cosmology commu- 
nity is now tasked with developing methods to analyze 
such data, and to extract as much information from 
various astronomical phenomena. These methods need 
to achieve unprecedented precision if the promise (and 
potential statistical power) of future surveys are to be 
fully exploited. 

Developing image analysis tools requires realis- 
tic mock data, containing as many instrumental ef- 
fects as possible, plus a known, underlying cosmo- 
logical signal, against which measurements can be 



judged. To generate galaxy sh apes in simulated 



ground-based images, Skymaker (lErben et all [2001) 



uses a simple physical model of concentric isophotes 
with a de Vaucouleurs profile for elliptical galaxies, 
and an additional, exponential component for spirals. 
By varying the model parameters, one can generate 
an unlimited number of unique simulated galaxies. 
However, deep field images from space-based tele- 
scopes contain galaxies with features more complex 
than these smooth analytical models can reproduce. 
Here we present t he full simage pipe l ine (as grad- 
ually d eve loped in IMassev et al. 20041 iMassey et alJ 
l2007bl and lFerry et al.l l2008). which empirically mim- 
ics the complex morphologies of galaxies seen in real 
data, such as t he H ubble Ultra Deep Field (UDF: 
Bec kwith et al. 2006) or Hubble Space Tel escope 



(HST) COSMOS survey (IScoville et al.ll2007l) . The 
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galaxy morphologi es are captured via a shapelet de 
composition (Refregier et al. 2003; Refregier & Bacon 



2003 



2005 



Bernstein & Jarvis 2002t IMassev & Refregier 
Mas sey et al. 2007al) . which also makes it easy 



to introduce a specified weak gravitational lensing sig- 
nal, useful to hone shear measurement methods. An 
example of shapelet galaxy im age simulation includes 
Sky lens dMeneghetti et al.ll2008h . 
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The simage code is written in the Interactive 
Data Language (IDL) and can be downloaded from 
|http : / / www . astro . caltech . edu/ ~r jm/ shapelet 
It also requires the core shapelets package, avail- 
able from the same location, and certain routines 
from the NASA Goddard Space Flight Center (GSFC) 
IDL astronomy user's library. Those required rou- 
tines are bundled in the simage download, but up- 
dated versions may be periodically available from 
|http : //idlastro . gsf c . nasa . gov/| In ad- 
dition, exploiting the full potential of some sec- 
tions of simage require SExtractor, available from 
|http : //astromatic . iap . f r/software/sextract 



Original UDF image Shapelet reconstruction 



Fig. 1. — An example of a spiral galaxy (from 
the Hubble Deep Field) modeled using shapelets 
(Mas sev & Refregierll2005l ). 



This paper highlights the full capabilities of the 
code, its general structure, and the a number of pos- 
sible uses. In section 2, we briefly detail previous ap- 
plications and associated results, plus the strengths and 
weaknesses of the utilized shapelet formalism. In sec- 
tion 3, we describe the structure of the software pack- 
age, and discuss the main modules of the software and 
their uses. In section 4, we exploit the software to in- 
vestigate telescope survey depth. We conclude in sec- 
tion 5. 

2. Capabilities and applications of 'simage' 

At the heart of simage is the ability of the shapelets 
method to efficiently and flexibly reconstruct com- 
plex galaxy morphologies. Shapelets are a complete, 
orthogonal set of basis functions, and a weighted lin- 
ear com bination of these can represent any localised 



image (Refregier et al. Il2003t iBernstein & Jarvisl[2002l 
(BJ02); Massey & Refr egier! \200$) . This is analo- 



gous to a Fourier transform, where weighted com- 
binations of sines and cosines can be used to recon- 
struct non-localised images. The mathematical prop- 
erties of the shapelet basis set make it particularly 
convenient for astronomical image processing, includ- 
ing quick convolutions with Point Spread Functions 
(PSF), pixillation, and operations such as transla- 
tions, rotations, magnifications and shears, that can 
be used to add a known signal into a simulated image 



dRefregier & Baconll2003t iMassev et al.ll2007al) . The 
ability to represent shears in a simple manner, allows 
for the inclusion of distorting shears produced by both 
the telescope's optical systems and weak gravitational 
lensing within galaxy clusters. While a Gaussian- 
based shapelet representation of a galaxy with a par- 
ticularly strong central peak, or extended wings, is not 
ideal due to the difficulty reproducing such features 



with Gaussian forms, a number of tests to reproduce 
exponential radial profiles with shapelet basis func- 
tions have show good reproduction of the majority 



of ga laxy shapes with very little bias ( IMassev et al 
20073). 



The simulated images are capable of being pro- 
duced in the B, V, i and z bands, since they are based on 
a galaxy morphology catalog pre-constructed from the 
Hubble UDF. Hence, the available passbands allowed 
by the F485W, F606W, F7775 W, and F850LP filters of 
the HST 

For a more in depth description of the method, in- 
cluding a general mathematical introduction to its ap- 
plication to image simulation, the reader is directed 
to the aforementioned references and Appendix A. In 
this section, we shall highlight the capabilities such a 
shapelet formalism provides, along with associated ap- 
plications and results. 

2.1. Prior applications: Shear studies, data codec, 
and mission development 

Previous applications of the shapelets image simu- 
lation, in it's simage incarnation as presented here, or 
in other forms, have ranged from direct image creation 
for community shear-analysis studies, data compres- 
sion investigations, and telescope/instrumental devel- 
opment. We will briefly detail these below. 

The Shear TEsting Program (STEP) was a collab- 
orative project which aimed to improve the accuracy 
and reliability of all weak gravitational lensing mea- 
surements in preparation for the next generation of 
wide-field surveys. STEP was launched in order to test 
and improve the accuracy and reliability of all these 
methods through the rigorous testing of shear mea- 
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Fig. 2. — Left: Plot of n e ff galaxy counts vs. exposure time at a pixel scale of 0.1 arcseconds. The n e ff quantity in 
this context is the number of galaxies within a given survey area that meet a particular desired size and signal-to-noise 
ratio. Right: Plot of n e f f galaxy counts vs. pixel scale at an exposure time of 400 seconds. These particular plots were 
generated from images created in the /-band (band 2 in simage) for a mock-up telescope 1.2m diameter mirror over a 
6 arcminute 2 area of the sky - see Table 1 . 



surement pipelines, the exchange of data and the shar- 
ing of technical and theor etical knowledge wit hin the 
weak lensing community (IMassev et al . 2007b). Here, 
the shapelets based simage code was used to create 
image data with incorporated shear fields for analy- 
sis in the testing program. Subroutines of the code 
were also used in the a nalysis of the resultant data. 
(ISchrabback et al. 2009) used images with simulated 
shear produced by this software to verify and vali- 
date their weak lensing code, which they then used 
to provide the first direct detection of dark energy us- 
ing weak lensing tomography. Opening up the is- 
sue to wider participation, particularly computer scien- 
tists, the GRav itational lEnsing Accuracy Testing 2008 
(GREAT08 : iBridle et al.l2008TlBridle et ai]20 09). and 
again allowed for the application of simage's image 
simulation capabilities. 

In Vanderveld et al. ( 2010l) . simage was utilised to 
test compression-decompression (codec) algorithms 
and methods for future visible survey telescopes. This 
is of vital importance given the vast quantities of data 
both space- and ground-based survey telescopes are 
predicted to produce in the comi ng years ( 10s of 
petabytes; e.g. LSST; llvezic l2007h . The roll of sim- 
age here was to create batches of simulated image data 
that recreated sizes, morphology, fluxes, and shapes 
that accurately mimic those we might expect from a 



visible survey telescope in both orbit and on Earth. 

Specific mission development has also been an 
application of the simage routines. In particular, 
the simulation pipeline has been a critical opti- 
mization tool in the optical design and observation 
strategy for ESA's co smic visions candidate Euclid 
( Refregier et al. 2010l) . Nasa/DOE Joi nt Dark Energy 
Mission (IDEM) concepts (e.g. SNAP; |jelinskvl 2006; 
High et al.l 120071) . and the weak lensing balloon mis- 
sions. 

2.2. Example application 

As a demonstration of one of simage's potential 
uses, we present an example investigation intended to 
explore telescope survey depth, specifically the rela- 
tion between the effective number of galaxies observed 
in a survey sample, n e ff, vs. pixel scale and expo- 
sure time for a mock-up space telescope design (the 
n e f f quantity in this context is the number of galaxies 
within a given survey area that meet a particular de- 
sired size and signal-to-noise ratio). A list of possible 
input telescope and instrument parameters are shown 
in Table 1 . We generate images according to these in- 
puts using a varying the pixel scale between 0.05 and 
0.40 arcseconds/pixel (with constant exposure = 400s), 
and the exposure time between 100-800 seconds (with 
constant pixel scale = 0.1'Vpixel). An example of the 
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Table 1 : Demonstration inputs within telescope. param 
for a mock-up telescope design 



• num_counts_frac.pro Calculates the galaxy 
magnitude distributions normalized t o the COS- 



MOS survey data at mid magnitudes (iLeauthaud et al 



Parameter file input: 


Value: 


throughput .ratio 


[0.0, 0.0, 1.2, 0.0, 0.0, 0.0] 


pixel.scale 


[0.0, 0.0, 0.1, 0.0, 0.0, 0.0] arcsec./pixel 


reacLnoise 


5 electrons 


collecting.area 


1.13m 2 


bancLbegin 


2 


bancLend 


2 


exposure.time 


100 seconds 


area 


6.0 square arcminutes 


n_star 


1000 


n_gal 


Default 


ee50 


0. 15 arcseconds 



output as created from a set of images simulated from 
simage is shown in Fig. [2] 

3. The simulation pipeline 

This section will detail the critical internal routines 
of simage, and describe the flow of the simulated im- 
age creation. 

3.1. Module overview 

The routine simage.pro is the primary program; it 
utilizes various routines within the pipeline to manu- 
facture a simulated image. Keywords listed in Table 2 
can be used to specify the desired telescope and survey 
characteristics. By default, the image will be produced 
in the B, V, i and z bands, based on a galaxy morphol- 
ogy catalog pre-constructed from the Hubble UDF. 
More permanent changes to the telescope and survey 
characteristics can be fixed in the telescope.param 
file. 

Figure[3]details the pipeline's main processes in the 
form of a flow chart. We note from the chart that 
there are three main stages to the pipeline; the multi- 
wavelength catalog generation, the repopulation of the 
catalog images into a field/resolution governed by the 
desired telescope parameters, and finally the addition 
of the various noise components. Other important rou- 
tines in the pipeline are, 

• simage_make_shapelet_object.pro: Generates 
a pixellated image of one simulated galaxy from 
a given a set of shapelet coefficients. 

• simage_make_analytic_object.pro: Generates 
an object for the image simulations, using an an- 
alytic profile. The size, magnitude and elliptic - 
ity are drawn from a real UDF galaxy template. 



20071) . and to a compilati on of other surveys at 



low and high magnitudes dMetcalfe et al. 2001 
see §3.2.2 for discussion.) 

• get_telescope_psf*: Reads in the desired PSF 
fits file before converting it into shapelet space 
(* 'telescope' represents a variety of telescope 
or survey names included in the pipeline e.g. 
getjudfjpsf, etc.). The PSF is an array of odd 
dimensions and be in logarithmic units, by con- 
vention. 

The desired size of the PSF is quantified in two 
ways: via the energy that encloses 50% of the PSF 
energy (EE50), or the Full- Width at Half Maximum 
(FWHM). The EE50 term is more commonly used 
in optical engineering, whereas FWHM is more com- 
monly used in science applications. Both are avail- 
able as user inputs and both have particular advantages 
when creating simulated images. For example, when 
comparing the EE50 to that of the FWHM, the former 
is more affected by the tail's profile. For this reason, 
EE50 is a better measure of how compact something is 
if you only have one number to describe an object. As 
a comparison, a Gaussian with a FWHM of 1 .77 pix- 
els, would have an EE50 of 0.885 pixels, while a typ- 
ical PSF with the same FWHM of 1 .77 pixels, would 
have an EE50 of 1.21 pixels. 

3.1.1. Input catalogue generation 

(|2008h . the first task 



As described in Ferry et al. 



is to generate a catalog of galaxy morphologies from 
real data. Note that galaxy morphologies are already 
provided from the UDF and, until better data become 
available, this time-consuming section of the pipeline 
need not be re-run. However, If desired it is possible 
to regenerate the UDF catalogue, or indeed to generate 
additional catalogues. 

The simulated images a re based on UDF data, with 
photometric redshifts from lCoe et al. ( 20061) . The pro- 
cess of getting from real data to simulated images, us- 
ing modules included in the simulation pipeline soft- 
ware package, is as follows. First, objects in the real 
data are detected and cataloged using the SExtractor 
routine on the image using a specified configuration 
file, config.sex an example of which can be found in 
the analysis directory of the simage package. The cat- 
aloged objects are then decomposed into shapelets by 
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Table 2: Descriptions of the user inputs for tele- 
scope.param. These parameter inputs govern tele- 
scope design and in turn the resultant output im- 
ages. Some parameter inputs are degenerate, e.g. if 
a throughput j-atio is entered, the filter Jiles path is ig- 
nored. 



Parameter file input: 


Description: 


throughput_ratio 


Total system throughputs relative to UDF 


sky.level 


Skylevel for each band (counts/s/arcsec 2 ) 


zeropoint 


Chosen zeropoint for the images in each band 


pixel.scale 


The instrument pixel scale in arc second/pixel 


psf.type 


Selects which PSF to use 


psf.path 


Path to the user's chosen PSF .fits file 


read_noise 


CCD read noise in number of electrons 


background 


= background subtracted, 1 = added 


collect ing_area 


The mirror collecting area in m 2 


bandjoegin 


The band on which to start the simulations 


band.end 


The band on which to end the simulations 


exposure-time 


Exposure time in seconds 


area 


The area on the sky to simulate in sq. arcmins 


random_seed 


A random seed for all random selections 


gamma 


The user specified weak lensing shear 


out put _f ile.pref 


Selection of output image file names 


n.star 


Number of field stars to be added 


n_gal 


Number of field galaxies 


f ilter.f iles 


Path to user's transition filter files 




The half light radius of the PSF 


fwhm 


The full-width at half-maximum of the PSF 



Table 3: Meaning of flags output from shex.pro detail- 
ing to the user how well a given object was modeled 
with shapelets. A value implies a successful decom- 
position, while a 10 signals failure. A FWHM of 
implies a profile fit to the object failed. 



Flag: 


Status: 





OK 


1 


Nearby object 


2 


Severe overlapping with nearby object 


3 


Object is near a saturated pixel 


4 


Object is near a masked region 


5 


Object is near the edge of the image 


6 


Object is itself masked out 


7 


Object has FWHM 


8 


Too few background pixels around object 


9 


Object entirely overlapped by neighbors 


10 


Routine sexcat2pstamp crashed 



running shex.pro. This program takes the real im- 
age, the SExtractor catalog, and the desired n max as 
inputs, and outputs a catalog of shapelet coefficients 
for the objects in the SExtractor catalog. We chose 
n max = 20 for the optical band - which is sufficient to 
model the HST PSF. shex.pro also uses the shapelets 
focus suite of routines to optimize the n max , /?, and 
centroid parameters. The program also outputs flags, 
from to 10, to tell the user how well a given object 
was modeled with shapelets, being good and 10 sig- 
naling failure. A list of the flags and their correspond- 
ing criteria is shown in Table 3. There is also an op- 
tion in shex.pro to remove a specified constant PSF, 
which can be modeled from the stars in the image. 
The PSF removal from the UDF catalog is a three step 
process. Firstly, the HST PSF is modeled by select- 
ing stars in the UDF images and decomposing them 
into a sum of shapelet basis functions. These stars 
are the best representation of the PSF contained within 
the image. Secondly, the galaxy objects of an image, 
once also decomposed into shapelet space, have the 
star/PSF shapelets subtracted, thus leaving a catalog 
of galaxies with the HST PSF removed. The stars in 
the newly created catalog are then discarded, since any 
PSF model that will be introduced will be formed from 
new simulated stellar images. 

Galaxies then have to be cross-matched across 
bands of data. This is done using srcor.pro, in the 
IDL Astronomy Library, to a tolerance of 1 arcsecond. 
Some galaxies will appear in all bands and some will 
appear in a subset, but each galaxy is given an ID num- 
ber and then a master catalog is created that contains 
the information for each unique galaxy ID, including 
its position, redshift, and shapelet coefficients in each 
band. For the UDF, this catalog is stored as a struc- 
ture called shapecat_total_trim.sav, which should be 
located in the specified data directory alongside the 
necessary psf folder. The simage program then ran- 
domly draws galaxies (in the form of their decom- 
posed shapelet coefficients) from this catalog when 
producing simulated images. 

3.1.2. Noiseless image simulation 

The pipeline would then proceed as follows: The 
simulated image is scaled according to the instrument 
and filter throughput. These are defined by parameters 
throughput j-atio and filter -files. The first specifies the 
throughput ratio for each band compared to that of the 
HST-Advanced Camera for Surveys (ACS), should it 
already known by the user. The second option allows 
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you to calculate this throughput j-atio array by speci- 
fying an arbitrary filter curve. This filter Jiles lists the 
total desired instrumental throughput for each band, in 
the format of wavelength (A) vs. throughput, where 
1.0 represents 100% transmission. The filter curve is 
then integrated and compared against the HST/ACS 
filter curve to generate the correctly normalized array 
throughput j-atio. In this way, it is possible to use UDF 
images but normalized for the throughput of any given 
telescope design. However, note that while the inte- 
gral of the throughput is considered, its shape is not. 
For example, the shape of galaxies is therefore not 
changed by a transmission curve that peaks redward 
of the HST/ACS filters and should therefore enhance 
the bulges of galaxies. The pipeline then reads in the 
specified PSF file, which remains constant throughout 
the simulation. 

The simage routine will then pass to simage .assemble 
Jmage which reads in the UDF shapelet catalog and 
populates the image with either n_gal number of galax- 
ies, or a pre-defined default that is calibrated to ex- 
isting observationa l data (the HST Cosmos survey; 
Scovil le et al. I l2007h . At this point any specified weak 
lensing shear, represented by the 2D array [71,72], 
is added to each galaxy, i.e a constant value across 
the field (see §3.4). The pixel scale of the shapelet 
catalog galaxies are adjusted to the specified value at 
this stage. A similar process is performed to pop- 
ulate the image with field stars, each of which are 
represented by the PSF. Here the subroutine sim- 
age ^tar jnagnitude distribution generates a random 
flux level for stars in the image. The stars follow the 
stellar luminosity fu nction measu red at the galactic 
poles and tabulated in lAllenld2000l) . 

We note that the pipeline also has the ability to cre- 
ate a simulated image containing objects with analytic 
(e.g. de Vaucouleurs) profiles with the same size, mag- 
nitude and ellipticity distributions as the UDF shapelet 
catalog. The shapelet coefficients are not used for any 
purpose other than the determination of these distribu- 
tions. Hence the pipeline can also use a non-shapelet 
method to create simulations, which can themselves be 
used to test shapelet based shape-measurement tech- 
niques if desired. 

The output noiseless image is written in units of 
photons or counts per second, along with a mock SEx- 
tractor output file with the all the objects' known 
input positions, sizes, magnitudes, ellipticities and 
star/galaxy classifications. All files are output to the 
Data directory. 



(A) Input: Real images e.g. Hubble Deep Field (HDF), Ultra Deep Field (UDF) 



shex.pro 



Galaxy selection: SExtratractor, on combined multi-wavelength stack 



Star selection 
Shapelet model of PSF 
Deconvolution 
Shapelet model of galaxy 



Catalog matching: simage_generate_shapecat .pro 



Output: 

Multi wavelength galaxy catalogs containing morphology, magnitude information etc. 



(B) Input: Multiwavelength catalogues, PSF model, simulated survey area 



simage_assemble_image . pro 



Decide number counts of stars 
Decide locations of stars on sky 



Determine locations of stars 
within image 



Pixellate stars (stars are the PSF) 
Add stars to image 



Decide number counts of galaxies 
via relation N(m) = B*10l A * m > 



Simple galaxy morphology: 



Decide locations of galaxies on sky 
Pixellate galaxies 
Add galaxies to image 



Save input parameters in FITS header, write to disc 



Output: 
Noise-free simulated image 



(C) Input: Noise-free image, exposure time, collecting area, throughput etc 



simage_add_noise . pro 



Read in noise-free image 



Add photon counting shot noise 



Add sky background level and Gaussian shot noise 



Add CCD detector effects, including CTE, dark current, etc. 



Add read noise 



Update FITS header, write to disc 



Output: 
Final simulated image 



Fig. 3. — Flow chart of the simulation pipeline's key 
processes beginning with the initial input UDF images 
to the final output .fits images. 



6 



E 



o 
o 

o 



200 



150 



100 



50 - 








0.0 0.5 1.0 1.5 2.0 

Galaxy Size - FWHM cut [arcseconds] 

Fig. 4. — A comparison of galaxy number density vs. 
galaxy size for both real HST-CO SMOS survey data 
(solid line; iLeauthaud et al. 2007) and simage simu- 
lated COSMOS data (dashed line). The galaxy size is 
measured via the SExtractor FWHM_ IMAGE. The x- 
axis scale is the FWHM cut size - i.e. galaxies that 
have a given size and above. Simulated data used 
COSMOS survey parameters of 0.03 arcseconds/pixel 
and an exposure time of 2000 seconds. 



Although the galaxy catalogues are created from 
the UDF, the magnitude distribution of the resultant 
simulated images is normalized to a variety of exist- 
ing galaxy surveys (described below) rather than just 
drawing randomly from the UDF galaxy magnitudes. 
Since the decomposed shapelet objects have color rep- 
resentative of real galaxies, this approach can apply to 
all bands available to the simulation pipeline (i.e. B, V, 
i and z). The routine num_counts_frac.pro utilizes the 
following number (N) counts relation; 



N = B* 10 (A * m) 



(1) 



where A and B are normalization factors, and m is the 
galaxy magnitude. The form of this expression and the 
values of the normalization factors A and B are taken 
from existing COS MOS survey analysis b etween mag- 
nitudes 21 and 26 dLeauthaud et alj|2007l) . COSMOS 
is the widest HST survey and as such is least affected 
by cosmic variance/sample size. It is also the closest 
data to the future space surveys that simage aims to 
model. Below magnitude 21 and above magnitude 26 



the the number counts are fit to var ious survey sources 
as compiled in lMetcalfe et alJ (1200 II) . The form is con- 
tinuous at these transition points and results in a real- 
istic number count relation for the simulated images at 
low, mid, and high magnitudes. 

In figure[4]we display a comparison of galaxy sizes 
in the HST-COSMOS survey and those from an simage 
simulated COSMOS field. We see that the normal- 
ization of simage to the COSMOS counts discussed 
above allows the pipeline to obtain a very similar sizes 
distribution. We note that there are differences at the 
lower size limit due to simage reproducing fewer very 
small galaxies in the simulation. This is because the 
pipeline does not decompose the smallest galaxies into 
shapelets since some of them are noise (or noisy) and 
as such do not end up in the shapelets catalog. 

3.1.3. Addition of noise 

The routine simage_add_noise.pro reads in the .fits 
format noise-free image, and writes out a noisy im- 
age, plus an inverse variance weight map. This is very 
fast to run, and is intentionally kept separate from the 
previous sections because a common task is to inves- 
tigate the effect of changing the survey exposure time. 
In this case, the same noise-free image can be used, 
and simage_add_noise.pro run multiple times in iso- 
lation, with different input parameters. Specific noise 
features and detector effects that can be added to an 
image post-creation, e.g. dark current, are included 
within the simage_add_noise.pro routine itself and are 
activated by selecting the corresponding flags when the 
routine is initially called. 

Noise is then added to the image. For convenience, 
the noise model is calculated in two separate compo- 
nents: shot noise on astrophysical sources, and on the 
sky background. In both cases, a random distribution 
of uncorrelated pixel values is drawn from a Gaussian 
with width equal to the square root of the counts in 
each pixel. The sky background level is estimated by 
default from that in the UDF, but can also be specified 
in telescope.param, in units of counts/second/arcsec 2 . 
By default, the constant sky background level is then 
subtracted, although this behavior can be turned off via 
the BACKGROUND keyword. 

Additional options include the ability to add read 
noise, to correlate the background noise to cheaply 
mimic the effects of DRIZZLE, to truncate satu- 
rated pixels, and to truncate at zero any non-physical, 
slightly negative pixel values (arising from noise or 
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modeling problems in the original UDF). None of 
these are enabled by default. To mimic the effects 
of DRIZZLE, the image is convolved with a kernel 
similar to the drizzle drop kernel and the original pixel 
square. This has the effect of correlating adjacent 
pixels in the image, and is most noticeable in the back- 
ground noise in blank areas well away from sources. 

4. Summary 

We have introduced an IDL simulation pipeline, 
simage, that creates mock deep field survey images by 
drawing from a catalog of Hubble Ultra Deep Field 
galaxies. Each galaxy in the catalog is decomposed 
into a set of analytic shapelet basis functions which can 
completely describe their morphological properties in 
a simple manner. We have shown how this catalog can 
then be used to populate a field of any given size and 
resolution depending on the user's requirements. The 
pipeline allows the user complete control over parame- 
ters such as exposure time, PSF type, mirror size, pixel 
scale, field star density, and noise, and simulates fields 
in the B, V, i and z bands. 

The code also has the ability to introduce a weak 
lensing signal into the data, allowing the output to 
be used for studies into weak lensing reconstruc- 
tion analysis. It is envisioned that the code will 
be used as a tool for research, instrumental de- 
velopment, and teaching. It is available to down- 
load as a self contained package of IDL modules at 
www.astro.caltech.edu/ ^rjm/shapelets. 
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Appendix A: Shapelets formalism 

Shapelets come in two flavors: Cartesian shapelets 
are separable in x and y, and polar shapelets in r and 
9. There is a one-to-one mapping between the two, so 



without loss of generality, we shall adopt whichever 
has the more convenient symmetries for the task at 
hand. The polar shapelet basis functions 
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where L^{x) are the Laguerre polynomials, have an 
overall scale size (3 and are parameterized by two in- 
tegers, n and m, which are the number of oscilla- 
tions in the radial and tangential directions. The ba- 
sis functions are calculated using shapelets_chi.pro. 
Using shapelets_decomp.pro, a galaxy (or star) im- 
age 7(r, 9) can then be decomposed into (complex) 
"shapelet coefficients" f n . m 



fn, m = J J I(r,6) X n,m(r,d;p)rdrd9, (3) 

so that the (wholly real) image can be reconstructed, 
using shapelets_recomp.pro as 



I(r,9) = 



oo n 

E E 

n— m— — n 



fn,mXn,m(r, 0; (3) . (4) 



In practice, it is necessary to truncate the expansion 
at some maximum value of n. Figure Q] shows an ex- 
ample galaxy image and its reconstructed counterpart 
using shapelets up to order n max — 20. It can be seen 
that the model easily captures the major features of the 
original galaxy. 

In shapelet representation, convolution between 
two images (such as a galaxy and a telescope's Point 
Spread Function) is simply a matr ix multiplicatio n of 
their /„, m coefficient arrays (Refregier et al. 2003). In 
the code, this is implemented via shapelets_convolve.pro. 
It is also possible to perform a deconvolution by in- 
verting the PSF matrix; this is incorporated within 
shapelets_decomp.pro. 

While previous operations were performed in po- 
lar shapelets since the functions are separable in r 
and 9 (rendering many operations more intuitive), 
pixillation can be performed most easily by switch- 
ing to Cartesian shapelets, then switching back. A 
closed form for the integrals of Cartesian shapelet ba- 
sis functions over rectangular pixels is given in §4.3 of 
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Massev & Refregierl ( 12005b . and is enabled by default 
in shapelets.chi.pro. 

We shall use several transformations of the galaxy 
images, first to randomize their appearance in the final 
simulation, and then to impose a gravitational lensing 
signal. In shapelet space, galaxies can be easily rotated 
by adjusting the phase of their (complex) coefficients 
f n . m , or reflected in the x-axis by taking their complex 
conjugates. A weak gravitational lensing shear signal 
7 can be applied to first order by mixing adjacent co- 
efficients according to the mixing matrix 
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y/(n-m + 2)(n - m + 4) f n +2,m-2 

(n - m){n - m - 2) f n -2, m +2 
y/(n + m + 2)(n + m + 4) /, 



n+2,m+2 



as described in §2.3 of Massev et al. 



d2007ab . which 
also provides similar operations for flexion. In the 
above, 1 corresponds to the identity operator, while 
S corresponds the to shear operator. Routines to im- 
plement such operations in practice are located in the 
shapelets/operations/ subdirectory. 

The above prescription for shear is only accurate 
to order 7. This will be insufficiently accurate for 
very high precision work, or if the gravitational lens- 
ing signal is particularly large. A new innovation for 
simage is that this transformation can now be general- 
ized to include higher order r ) 2 , 7 s , etc. terms. This is 
achie ved mathematically by e xponentiating the opera- 
tion (Bernstein & Jarvisll2002h . For a practical imple- 
mentation to fourth order, note that the first four terms 
in an exponential expansion 



—2 ^3 ^4 
f 9 S S S 3 

1 + 5+ ^! + 3! + 4! =8 



1 + S {1 + S) 2 (1 + S) 4 



24 



(6) 



where S here is the shear operator but could equally be 
replaced by any other. To perform this on a shapelet 
model we simply need to apply the linear mapping © 
four times, recording the new coefficients f n m at each 
stage, and add them to the original coefficients in the 



, | . . A. This behavior is controlled 
via the ORDER keyword in shapelets_shear.pro, and 
is set to 4 by default. 

Note that, as discussed in BJ02, there is a some- 
what arbitrary choice for these higher order terms, 
which can be changed depending on the required def- 
inition of shear. The expansion above changes an in- 
trinsically circular source into an ellipse with major 
and minor axes a and b via a distortion 5 = (a 2 — 
b 2 )/(a 2 + b 2 ). A "conformal shear", v = arctanh(5), 
produces a slightly different ratio of major and minor 
axes, but can be achieved by simply adjusting the in- 
put shear. A fourth-order implementation of a confor- 
mal shear in shapelet space perfectly matches the real- 
space transformation of highly oversampled images 
within a computer's numerical precision (B. Rowe, 
priv. comm. 2008). It is therefore not just faster, but 
should be accurate within 1 % for shears up to 7 ps 0.47 
dMassev & Goldbergll2008l) . For a typical cosmologi- 
cal gravitational lensing signal of a few percent, ap- 
plying only a first order shapelet-based shear yields 
pixel values in the final image that are incorrect at a 
level of approximately one part in 10~ 3 , and chang- 
ing from 8 to v yields differences of around one part 
in 10 -5 . Because of this, the simage pipeline routines 
such simage_make_analytic_object.pro, uses the con- 
formal shear, v, when called to include a shear signal. 
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